Bray‐Curtis (AFD) differentiation in molecular ecology: Forecasting, an adjustment ( A A), and comparative performance in selection detection

Abstract Geographic genetic differentiation measures are used for purposes such as assessing genetic diversity and connectivity, and searching for signals of selection. Confirmation by unrelated measures can minimize false positives. A popular differentiation measure, Bray‐Curtis, has been used increasingly in molecular ecology, renamed AFD (hereafter called BCAFD). Critically, BCAFD is expected to be partially independent of the commonly used Hill “Q‐profile” measures. BCAFD needs scrutiny for potential biases, by examining limits on its value, and comparing simulations against expectations. BCAFD has two dependencies on within‐population (alpha) variation, undesirable for a between‐population (beta) measure. The first dependency is derived from similarity to GST and FST. The second dependency is that BCAFD cannot be larger than the highest allele proportion in either location (alpha variation), which can be overcome by data‐filtering or by a modified statistic A A or “Adjusted AFD”. The first dependency does not forestall applications such as assessing connectivity or selection, if we know the measure's null behavior under selective neutrality with specified conditions—which is shown in this article for A A, for equilibrium, and nonequilibrium, for the commonly used data type of single‐nucleotide‐polymorphisms (SNPs) in two locations. Thus, A A can be used in tandem with mathematically contrasting differentiation measures, with the aim of reducing false inferences. For detecting adaptive loci, the relative performance of A A and other measures was evaluated, showing that it is best to use two mathematically different measures simultaneously, and that A A is in one of the best such pairwise criteria. For any application, using A A, rather than BCAFD, avoids the counterintuitive limitation by maximum allele proportion within localities.


| INTRODUC TI ON
Comparisons of biodiversity between regions are important aspects of understanding both ecological and genetic systems. There are many geographic differentiation measures, used for purposes such as assessing genetic diversity and connectivity (Gruber et al., 2018;Guillot et al., 2005;Manni et al., 2004;Meirmans, 2020;Sherwin et al., 2017Sherwin et al., , 2021 and searching for signals of different selective regimes geographically, which is expected to have high false-positive rates (Bierne et al., 2013;Lotterhos & Whitlock, 2014;Narum & Hess, 2011;Schneider et al., 2021;Xiang-Yu et al., 2016). Because of the anticipated high false-positive rates, it is important to confirm findings using a wide range of mathematically unrelated measures. Often these measures are chosen from the Hill "Q-profile," which includes: counts or sharing of allelic types (Q = 0 measures); Shannon information and differentiation (Q = 1); and heterozygosity, nucleotide diversity, Gini-Simpson, G ST , F ST , Morisita-Horn, D EST (Q = 2) (Chao et al., 2014;Gaggiotti et al., 2018;Jost, 2008;Jost et al., 2010;Sherwin et al., 2017Sherwin et al., , 2021. However, despite their different sensitivity to some matters, such as rare and common alleles, the members of the Q-profile are all mathematically related (Sherwin et al., 2017(Sherwin et al., , 2021. Notably, one recent addition to the range of measures in molecular ecology is outside the Hill Q-profile: the Bray-Curtis index of dissimilarity, a method of assessing differentiation that is extremely popular in its original use, to assess differentiation between species assemblages (Bray & Curtis, 1957). During 2021 alone, Bray-Curtis was mentioned over 10,000 times in Google Scholar. Bray-Curtis can be expressed in a way that facilitates comparison with differentiation measures derived from Hill numbers; the mathematical equivalence to other formulations of Bray-Curtis is documented in (Chao & Chiu, 2016;Jost et al., 2010;Ricotta et al., 2021;Ricotta & Pavoine, 2022;Ricotta & Podani, 2017).
where a 1j and a 2j are the abundances (counts or frequencies 0 ≤ a ≤ ∞ ) in each of two locations (1,2), for variant j (1 ≤ j ≤ S) and S is the total number of species. This measure satisfies many of the requirements of a good measurement of differentiation between assemblages (Chao & Chiu, 2016;Magurran, 2004;Ricotta & Podani, 2017). This index is also used for analysis of operational taxonomic units in metagenomics (Peng et al., 2020).
Unification of ecological and genetic approaches is desirable, because of their interaction as parts of the same biological systems, and because of their underlying mathematical similarities (Rosindell et al., 2015;Sherwin, 2018), so it is good to see that a simplified version of Bray-Curtis has been proposed as a measure of differentiation in molecular ecology and evolution (Berner, 2019a(Berner, , 2019bShriver et al., 1997), echoing a similar measure in community ecology (Whittaker, 1975) (p 118). It was renamed "allele frequency difference" AFD, but I will call it BCAFD, in deference to its original proponents, and because it is a difference of proportions (0 ≤ p ≤ 1) rather than frequencies (0 ≤ a ≤ ∞). In the two-variant two-location case, Bray-Curtis simplifies to the unsigned difference between locations 1 and 2 of proportions of either of the two allelic variants (Berner, 2019a(Berner, , 2019b. where p 1 = a 1 ∕ a 1 + a 2 and q 1 = 1 − p 1 , and similarly for the other location p 2 and q 2 . When there are multiple alleles, it is suggested to use the sum of the absolute allele proportion differences divided by two ( [Berner, 2019a], Table S1), which actually is equivalent to the more general Equation (1). However, unless otherwise stated this article will deal with the biallelic case which is very common in current molecular ecology-SNPs or single-nucleotide polymorphisms.
For applications including selection detection and assessment of connectivity between locations, it is critical to know the measure's null behavior, that is, in the absence of selection ("neutrality"), with specified conditions such as population size, dispersal, and mutation (Bierne et al., 2013;Gruber et al., 2018;Guillot et al., 2005;Lotterhos & Whitlock, 2014;Manni et al., 2004;Meirmans, 2020;Narum & Hess, 2011;Schneider et al., 2021;Sherwin et al., 2017Sherwin et al., , 2021Xiang-Yu et al., 2016). Despite not belonging to the Hill Qprofile, BCAFD appears to have some mathematical relationship to two of the Hill measures: G ST and F ST (Appendix 1). Therefore, based on forecasts for those two measures, it will be shown that it is possible to develop forecasts for BCAFD for two-location, two-variant systems such as single-nucleotide polymorphisms (SNPs).
All diversity measures must be scrutinized for their particular properties (Leinster & Cobbold, 2012;Leinster, 2021;Sherwin et al., 2017;Sherwin et al., 2021). An important property of differentiation measures is independence between alpha (within location) variation, beta (between location) differentiation, and total (gamma)
Critically, G ST and F ST are well-known to have the serious limitation of being heavily influenced by within-location variation (alpha), something that is not desirable in a between-location (beta) differentiation measure. Although F ST was explicitly proposed as a measure of between-subgroup differentiation (Wright, 1943) and has been used for that extensively, unlike some other Hill-profile measures, F ST shows strong dependence on alpha within-locality diversity, as does the related measure G ST (Jost, 2008;Meirmans & Hedrick, 2010;Nei, 1977Nei, , 1973. Because of its relationship to G ST and F ST , it is likely that there will be dependency of BCAFD on alpha variation. Another dependency of BCAFD on alpha variation is that it is obvious from Equation (2) that BCAFD can never be larger than p max , the higher of the two allele proportions, p 1 and p 2 . In other words, if either p 1 or p 2 is zero, then the value of Bray-Curtis will be equal to the other, more abundant, proportion. Of course, the values p 1 and p 2 are within-location proportions of one of the two alleles-a within-population (alpha) measure. This is an extremely counterintuitive limitation on a between-location (beta) differentiation measure, and is expected to result in biased values. This might be particularly important when using the measure to search for loci that experience different directions of selection in different locations, because this difference of selective regime will obviously give a signal of large differentiation values between locations, relative to other neutral loci. As a result, the truncation of large values of BCAFD due to p max might be expected to reduce the ability to distinguish such adaptive loci from neutral loci.
The confounds with alpha variation due to relationship to G ST , and restriction by maximal allele proportion p max , require examination in this article; however another possible confound does not appear to be of concern. As well as the proportions of variants, a between location (beta) differentiation measure can be confounded by the number of variant types. This confound can be avoided by restriction to two-variant systems such as SNPs, as is done in this article. Also, it does not appear to be a problem for the multiallelic version of BCAFD (Equation (1), also [Berner, 2019a] Table S1). When there is maximal differentiation, that is, no alleles shared between locations, one expects to always get the maximal value for the genetic differentiation statistic. This in fact does happen. For example, if there are four alleles w, x, y, and z, with w and x in location 1, and the other two in location 2, so that p 1w = p 1x = p 2y = p 2z = 0.5, and other proportions are equal to zero, then the multiallelic statistic is equal to BCAFD = 1.0. Also, if location 1 only has allele w, and the other three alleles are in location 2, with p 1w = 1; p 2x = p 2y = p 2z = 1 3 , then the multiallelic statistic remains BCAFD = 1.0, as expected for the same situation of maximal differentiation (no shared alleles).
Irrespective of these confounds, it should be noted that the alpha-dependency of G ST ∕ F ST does not forestall all use of these measures, provided that we know their behavior under selective neutrality with specified conditions such as population size, dispersal, and mutation (Bierne et al., 2013;Gruber et al., 2018;Guillot et al., 2005;Lotterhos & Whitlock, 2014;Manni et al., 2004;Meirmans, 2020;Narum & Hess, 2011;Schneider et al., 2021;Sherwin et al., 2017Sherwin et al., , 2021Xiang-Yu et al., 2016). With this in mind, and responding to the increased use of BCAFD in molecular ecology described above, this paper carries out the following tasks: • It creates a modified version of BCAFD termed A A ("Adjusted AFD") that is corrected for the limitation by p max .
• Forecasts are made and tested, for Bray-Curtis (BCAFD) and A A, for selectively neutral biallelic SNPs-a very common data type at present-under various scenarios of population size, mutation, and dispersal. This will allow BCAFD, and especially A A, to be used for evaluating competing models of population connectivity, making projections for the future, or identifying outlier loci whose differentiation level departs from neutral expectations, and so are candidate adaptive loci.
• Simulations are performed to investigate how the A A correction for bias performs in detecting loci under directional selection, in comparison to competing measures, or in consort with those measures.

| MATERIAL S AND ME THODS
Forecasting equations for Bray-Curtis were developed for the common and simple case of a single neutral biallelic SNP locus, with two locations (1,2); the measure can be averaged over multiple loci, and can be applied to haploids, or to diploids in Hardy-Weinberg equilibrium (i.e., no population-wide correlation between the two alleles within diploid genotypes). When there are only two variants, the Bray-Curtis equation is: BCAFD = | | p 1 − p 2 | | (Berner, 2019a, 2019b) (Equation 2, above) where p 1 and p 2 are proportions of one of the two alleles at each location (q 1 = 1 − p 1 ; q 2 = 1 − p 2 ). (2) is a transform of two well-known differentiation measures (Halliburton, 2004;Wright, 1943):

The quantity in Equation
where 2 p is the variance of p between locations, H is the Hardy-Weinberg (Binomial) expected heterozygosity, for example, ; and p is the average p over the two identical in the two-allele, two location case ( [Halliburton, 2004] Box 9.5). Appendix 1 shows that Because of its close relationship to G ST or F ST , BCAFD forecasts can be based on well-known forecasts for those measures (Appendix 1). The expectation for diploid BCAFD at drift-dispersalmutation equilibrium is: (3 and A1.2) where m is symmetrical dispersal between the two locations (0 ≤ m ≤ 1); μ is the rate of mutation (0 ≤ μ ≤ 1); N is the effective population size at each location (identical); and 2 D is the second order Hill diversity, or effective number of alleles 2 D = 1 ∕ 1 − H T .

The equivalent equation for the haploid SNPs simulated in this
article is: The performance of these equations was assessed by simulation of biallelic neutral single-nucleotide polymorphisms (SNPs) in two haploid subpopulations, for a wide range of scenarios covering all possible combinations of three symmetric dispersal rates (m = 0.01, 0.03, 0.1) and three subpopulation effective sizes (N = 1000, 10,000, 100,000). Starting allele proportions in each subpopulation (p values) were randomized in each replicate. Simulations used the typical SNP mutation rate (μ = 10 −9 ), but essentially identical results were obtained with rates between μ = 10 −6 and 10 −12 . The simulation was programmed in MATLAB, and full details are in Appendix 2, and Dewar et al. (2011). There were 1000 replicate iterations of each scenario, which could also be considered as 1000 independently inherited loci (i.e., in linkage equilibrium). Each iteration was run for 200 generations, and each generation included stochastic binomial sampling of the parents' alleles to establish the allele proportions for the offspring, followed by symmetrical dispersal to create the parent populations for the next generation. Because the forecasts are for drift-dispersal-mutation equilibrium, it is important to know whether the simulations had reached equilibrium. The adequacy of the run-time of 200 generations was confirmed in three ways, detailed in Appendix 2: 200 generations was several times longer than the expected time to half-equilibrium values; inspection ensured an asymptote to a stable value for BCAFD; and the variance of BCAFD between-generations was much lower than variance between replicate iterations (typically one tenth or less). The performance of the simulation was checked by comparison with results of EASYPOP (Balloux, 2001) and with known predictions for G ST (see Appendix 2 for details).
To assess whether the expectation from Equation (6) was an adequate forecast of BCAFD, BCAFD was calculated at the final generation, then linear regression was used (in EXCEL). If the expectation from Equation (6) is accurate, it is expected that a regression of the simulated BCAFD against the expected BCAFD should have a slope of unity and an intercept of zero. Additionally, alpha-dependence was assessed, and possible corrections suggested, including an adjusted measure A A that has no limitation by p max .

In other investigations, I examined the relationship between
BCAFD and three other differentiation measures: G ST , D EST , and mutual information, I (Sherwin et al., 2017(Sherwin et al., , 2021. I also examined whether the forecasts could be made completely independent of within-location variation. Finally, I produced nonequilibrium forecasts, suitable for situations where there has been recent disturbance to connectivity, for example. Simulations were used to investigate the effect of the adjusted measure A A on detectability of loci under different directional selection in each population. These simulations were identical to the ones described above, with two alterations. First, the simulations were restricted to large population size and low dispersal (N = 100,000, m = 0.01). Second, selection was simulated each generation by, in one location, increasing the number of surviving progeny of one genotype by multiplying by a factor of 1 + s/2, and decreasing the same genotype by 1 − s/2 in the other location  (Jost, 2008); and mutual information I (Sherwin et al., 2017(Sherwin et al., , 2021. For each measure, I tallied the percentage of loci (out of 1000 simulated) that would be identified as outliers (i.e., potentially under selection) using the "univariate" criterion that their genetic differentiation values were in the top 1% of the 1000 loci simulated without selection in a parallel neutral simulation, separately for each one of the five differentiation measures. As well as those univariate criteria, the same analysis was repeated using a series of more restrictive "bivariate" criteria, that is, that for a locus in the selection simulation to be identified as an outlier, it was required to have differentiation in the top 1% of neutral loci for each of a pair of the differentiation measures

| RE SULTS
Trials of Equation (6) used the data from the haploid simulation program described above. Figure 1a  • there appears to be an oblique upper bound to the scatter of points from the 1000 replicates of each scenario; this will be discussed later.
• Despite the scatter of replicates, there is an extremely good regression of simulated BCAFD on predicted BCAFD (significance P was extremely low-assigned to zero by the program, see caption of Figure 1a). Note that the scatter is not unexpected given that the initial allele proportions were randomized.
• the intercept is extremely close to zero, as expected • however, the slope is slightly below the expected 45-degree line for perfect prediction, with a slope of 0.83, see caption of Figure 1a; the 95% confidence limits for the slope were 0.81 to 0.85, so that the limits did not include the expected unity.
In the introduction it was pointed out that the value of BCAFD is restricted by the maximum p value p max in either of the two locations, at the generation where BCAFD is calculated. This is a potential reason for the oblique upper bound for the observations in Figure 1a.
To Investigate this, the regression of simulated BCAFD on expected BCAFD was repeated on ten subsets of the 9000 datapoints, subdivided by the final value of p max , the maximum p in either of the two locations. Results in Table 1 show that the departure from a 1:1 slope is indeed due to the restriction by p max . The bottom two rows of this table are where there is the least constraint on simulated BCAFD values (0.8 ≤ p max ≤ 0.899 and 0.9 ≤ p max ≤ 1), and in these two cases the slope of the regression of simulated BCAFD on expected BCAFD is indeed unity as expected. The slope of this regression decreases linearly when it is more constrained, that is, with lower p max values (Table 1 and Figure 2).
There are two possible corrections for this dependency on maximum p value. First, the data could be filtered to only include loci with very high maximum p values (0.8 ≤ p max ≤ 1, Table 1, Figure 2), but of course this would greatly reduce the usable data. Second, because the regression in Figure 2 is very linear, one can correct the expectations for the effect seen in that figure, where (coefficient of simulated BCAFD on expected BCAFD) = 0.6152 + 0.3985 × p max , so that we create a modified version of BCAFD, called " A A" which is free of dependence upon p max : We then find that the forecasts are general for all values of p max , for haploid: or the same for diploid loci in Hardy-Weinberg equilibrium, replacing 4N with 8N: Figure 1b shows the plot of A A (i.e., BCAFD adjusted to compensate for limitation by p max ) plotted against the expectations from (Equation 8). This regression shows the expected slope of unity and intercept of zero, demonstrating that the simulation confirms the haploid prediction for A A in Equation (8), including for each individual scenario (Figure 1c).
There are nonlinear relationships between A A and three other differentiation measures: G ST , D EST , and mutual information, I, as was suggested by a previous investigation of BCAFD (Berner, 2019a(Berner, , 2019b ( Figure 3). This shows that A A provides information that is not linearly dependent on these other measures, which is important when using multiple measures for confirmation of results such as assessment of connectivity, and searches for loci potentially under selection.
As well as the equilibrium forecasts just described, it is important to have nonequilibrium forecasts for A A, which will often be relevant in many situations, including recently disturbed populations; nonequilibrium forecasts are shown in Equation (A1.11b).
It was also investigated whether the dependence of BCAFD on within location (alpha) variation could be fixed by basing the expectations for BCAFD not on G ST , but upon G" ST (Meirmans & Hedrick, 2010). Unlike G ST , G" ST is free of influence of withinpopulation variation. In Equation (A1.14), it can be seen that this new formulation of BCAFD is still heavily dependent upon heterozygosity H, including the within population (alpha) measures H 1 and H 2 .
With the false detection of selection held constant at 1%, the important matter is the performance value: what percentage of loci that are classified as outliers, due to their differentiation value surpassing the univariate or bivariate criterion, are actually under selection-the true positives (TP). For a wide range of selection strengths, Table 2 shows the performance values for each univariate criterion (a single differentiation measure), and each bivariate criterion (i.e., an outlier locus must surpass the cutoff value for two differentiation measures). Of course, with the strongest selection (s = 0.05), all criteria performed well, but with very weak selection (s = 0.001), there was poor performance. The right-hand column of Table 2 shows the performance averaged over all selection strengths, which had similar rankings for the performance of the criteria. The univariate criteria did not perform as well as the bivariates, with no overlap of mean performance ±1 × SE. Within the univariates, there was similar performance for all criteria, but when averaged over all selection strengths, the three best performers were A A, I, and G ST . Within the bivariate criteria again there was similar performance for all criteria.
Nevertheless consistently the three best performers were "

| DISCUSS ION
Science progresses by making forecasts under given conditions, then testing to see whether these conditions are confirmed by the data. Examples include assessing levels of dispersal by identifying whether neutral loci depart from expectations for isolation or panmixia, and testing for loci that may be responding to geographically variable selection, by identifying whether genetic differentiation is higher than neutral expectation ("outlier loci," (Bierne et al., 2013;Lotterhos & Whitlock, 2014;Narum & Hess, 2011;Schneider et al., 2021;Xiang-Yu et al., 2016)). Unfortunately, there are expected to be many false results in such molecular ecological methods (Bierne et al., 2013;Lotterhos & Whitlock, 2014;Narum & Hess, 2011;Schneider et al., 2021;Whitlock & McCauley, 1999;Xiang-Yu et al., 2016). Therefore, it is advisable to confirm conclusions by methods that are mathematically independent or at least partially independent. Figure 3 shows that A A =  (8) and (9). However, the alphadependence is not fatal; despite their alpha-dependence, G ST and F ST are frequently used in various ways, including assessing connectivity and searching for loci under geographically variable selection.
Moreover, under some conditions G ST and F ST have performance comparable or better than other measures (Schneider et al., 2021).
Nevertheless, like all such methods, there are expected to be many false-positives, so that corroboration with semi-independent assessments is needed (Bierne et al., 2013;Lotterhos & Whitlock, 2014;Narum & Hess, 2011;Schneider et al., 2021;Xiang-Yu et al., 2016), which is where A A might be used.
The neutral forecasts for A A can be used either to make biologicalinventories of differentiation between locations (or times), or to be compared to observations in order to assess biological processes that underlie all biology, and are the processes which some conservation initiatives aim to conserve (Anonymous, 1988). Processes to be investigated include population size, mutation, and dispersal in natural or managed systems, or searches for outlier loci that depart from neutral expectations, and are thus candidate adaptive loci, which of course are very important in evolution and conservation (Teixeira & Huber, 2021 Note: The 9000 data points from Figure 1a, sorted by p max in the final generation. In the first column, "Central p max = 0.05" identifies the points with 0 ≤ p max ≤ 0.099, etc. The remaining columns show the results of regression analysis of (Simulated-BCAFD) against (Predicted-BCAFD from Equation 6) for the subset of the datapoints identified in the left column. All regressions showed an intercept very close to zero, as expected. Large numbers of significant digits are retained in the slope coefficients because of their subsequent use in the analysis in Figure 2, where the coefficients are plotted against central p max values.
F I G U R E 1 (a) Comparison of simulation results with algebraic predictions for BCAFD; 9000 points from the 1000 replicates of each of nine neutral scenarios (effective size N = 1000, 10,000, 100,000, dispersal rate m = 0.01, 0.03, 0.10) and with regression equation (Simulated-BCAFD) = 0.83 × (Predicted-BCAFD from Equation 6) (significance P <<0, R 2 = .50, intercept negligibly different from zero: −7.6 × 10 −5 ). The black line is the regression line; the red line is the expected 1:1 relationship. (b) The same data again, using the correction for the limitation by maximum p, that is a plot of A A = |p 1 -p 2 |/(0.6152 + 0.3985 × p max ) against the expectation shown in Equations (6) and (8). In this case, the expected 45-degree plot is achieved exactly, with the expected slope of unity (slope coefficient = 1.00, 95% confidence limits 0.98 to 1.02, significance P <<0, R 2 = 0.50, intercept negligibly different from zero: 0.0004). The red line for 1:1 slope is exactly coincident with the regression line. (c) The nine scenarios from (b) plotted individually-comparison of simulation results with algebraic predictions, using A A, the correction for the limitation by maximum p. Each panel shows 1000 points from the 1000 replicates of one scenario, whose dispersal rate m and effective size N is shown in the panel's headline. The slopes of regression lines are shown on the panels, with 95% confidence intervals, which included unity in all except two marginal cases, and are therefore each concordant with the overall result shown in (b) and the relationship in Equation (8). In all cases, the intercept was negligibly different from zero, and P for significance was <10 −18 .   Table 2).
The differences in performance were small, but even small improvements are very important given that this commonly used approach can only identify outlier loci that are putatively under selection, then each of these "candidate" loci must be confirmed by separate extensive investigations, such as "evolve and resequence" experiments in one or more standard environmental conditions (Schlötterer et al., 2015).  (Jost, 2008;Magurran, 2004;Sherwin et al., 2017;Sherwin et al., 2021). If the primary purpose of assessing differentiation is for identification of loci under selection, another good measure to contrast with identifications by A A would be B GD , which can be used at any level of the Hillfamily "Q," and has a good sensitivity to selection, and is particularly appropriate for multi-SNP haplotypes, which are not considered in the current work (Schneider et al., 2021).
Of course, any use of theory relies upon adherence to assumptions, and this initial theory for A A has assumptions like any theory.
The equations for G ST , upon which the A A forecasts are based, rely on a number of assumptions (Neigel, 2002;Ochoa & Storey, 2021;Semenov et al., 2019;Whitlock & McCauley, 1999) and each of these needs to be investigated if it is proposed to apply Equations (8) or (9) to any particular case. is shown as squares, G ST as discs, I as triangles. All measures were from the same simulated dataset that was used in Figure 1.

F I G U R E 2
The effect of maximum p-value p max on the regression slope coefficient of (simulated BCAFD) on (expected BCAFD from Equation 6). This plot itself has a regression equation: (coefficient of simulated BCAFD on expected BCAFD) = 0.6152 + 0.3985 × p max , with R 2 = .90, and P = .000025. The values upon which the plot is based are taken from Table 1

± 4.40
Note: The table shows the number of loci (±SE) from selection simulations of 1000 loci, which were identified as being under selection by criteria based on differentiation values from neutral simulations of 1000 loci: either a "univariate" criterion of being in the top 1% of neutral values for one differentiation measure, or a "bivariate" criterion of being simultaneously in the top 1% for two differentiation measures. In each of columns 2-5, the top value in each cell is the number of loci identified as being under selection (true positive, TP), in the selection simulation with the known value of selection shown at the top of the column, out of the total of 1000 independent loci simulated. The second value in each cell is the number of loci identified as being under selection (False positive, FP), in the parallel neutral simulation; of course with univariate criteria and the cutoff being the top 1%, the FP value is always 10 (1% of 1000 loci). The third value in each cell is the "performance" value-the percentage of loci that are true positive, out of all loci identified as outliers by that criterion (TP & FP). The performance value shown is for the case where 1% of all loci were under that selective regime, and all other loci were neutral; the calculation is 100 × (TP × 0.01)/[(TP × 0.01) + (FP × 0.99)]. Of course, the proportions of neutral and selected loci would not be known beforehand in a study designed to detect loci under selection, but given that it is standardized to a constant univariate FD rate, the performance values can be used to compare the criteria. The right column shows the performance averaged over all four selection strengths. Within each of the univariate criteria and the bivariate criteria, the three criteria with the best average performance are bolded. Note that the rank order of performance values is similar for most selection strengths, except the weakest selection (s = 0.001). (Chao & Jost, 2015) (A. Chao pers. comm.); however, this correction method is also inapplicable to any measure based on a two-variant system such as SNPs. Finally, Figure 1b shows a wide scatter, but the regression analysis shows that if there are multiple independent replicates such as hundreds, or a thousand, neutral SNP loci in linkage equilibrium, the neutral forecast is very accurate. This number of statistically unlinked SNP loci is easily achievable with current methods for genotyping-by-sequencing (e.g., www.diver sitya rrays.com).
Irrespective of whether one wishes to use theoretical expectations, it is advisable to use A A rather than BCAFD, because the latter's dependence on p max limits its comparability to other studies, even within the same species, if the population pairs analyzed are in parts of the range that have different p max , due to a strong cline. if the speciation rate is negligible relative to the dispersal rate; this is of course the original use of Bray-Curtis (Bray & Curtis, 1957), which would require development of multivariant theory plus simulations tailored to species assemblages, including investigation of the wide scatter seen in Figure 1, for which species analyses could not be overcome by using hundreds or more replicate loci-instead, hundreds or more replicate pairs of communities would be needed, which is probably unattainable.
In conclusion: • The new A A measure (Equation 7) provides a semi-independent means for assessing connectivity, selection, etc. based on geographic genetic differentiation, that can be used in combination with other such measures to minimize errors such as false positives.
• The A A measure avoids counterintuitive truncation of high values of beta-differentiation by alpha within-population variation (p max ), • Avoiding this truncation means that that studies with different p max can now be compared realistically, either between species, or even within the same species, if the population pairs analyzed are in parts of the range that have different p max , due to a strong cline.
• Avoiding this truncation is especially important if the high values of differentiation are to be used to identify candidate adaptive loci, because the truncation would pull the truly high values in amongst the not-quite-so-high, leading to increased false negatives and positives.
• As predicted, the best performance at identifying outlier loci that are potentially under selection comes from using two geographic genetic differentiation measures simultaneously, to make bivariate criteria; the three best performers were "D EST & I," tied with The differences in performance are very important given that each of the identified "candidate" loci must be confirmed by separate extensive investigations • As well as simply presenting patterns in the data, if researchers consider that their system conforms to the assumptions herein, the neutral forecasts for A A can be used as a rigorous basis for investigations such as tests for selection and assessment of connectivity.
• There are equilibrium and nonequilibrium versions of the theory for A A (Equations 8, 9, A1.11b).
• Irrespective of whether the theory in this paper is used, BCAFD cannot be free of the limit of maximum within-population allele Warton, Jia Zhou, and an anonymous reviewer.

CO N FLI C T O F I NTE R E S T
There are no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
MATLAB program, and data, are on DRYAD at https://doi. i.e.,: In Equation (A1.11a), (A1.11b), for haploids, 2 N is replaced by N.

GAMMA?
Note that due to the dependence of G ST on (alpha) heterozygosity, BCAFD, and derived measures such as A A, are also expected to be dependent upon (alpha) heterozygosity, and such dependence is not a desirable property for a measure of between-locality (beta) differentiation. This is additional to the dependence of BCAFD on maximum allele proportion, for which a correction is applied in the main Although G ′′ ST is free of alpha-dependency, nevertheless this new formulation of A A is clearly still heavily dependent upon heterozygosity (H T , H 1 ,H 2 ), which are gamma and alpha measures, which would ideally not affect a (beta) differentiation measure. Additionally, using this formulation in Equation (A1.14) for BCAFD would considerably complicate the derivation of theoretical expectations comparable to Equations (A1.7) and (A1.8).

A PPE N D I X 2
The MATLAB simulation program This MATLAB program was modified from the one previously described (Dewar et al., 2011), to include calculation of the Bray-Curtis/ AFD Index BCAFD and A A (Equations 2, 7, A1.1), as well as the previ- .2), D EST , and I (mutual information).
The simulation deals with two biallelic haploid subpopulations, for scenarios with every possible combination of levels of three symmetric dispersal rates (m = 0.01, 0.03, 0.1) and three effective subpopulation sizes (N = 1000, 10, 000, 100, 000), giving a total of nine scenarios. Mutation rate between the two alleles per SNP locus per generation was held at μ = 10 −9 , a likely rate for SNP alleles, but preliminary simulations showed virtually identical results with the much higher and lower rates of μ = 10 −6 and 10 −12 . Indeed, inspection of Equations (A1.7) and A1.8 show that mutation rate will be largely irrelevant unless the locations are completely isolated (m = 0 exactly).
Starting allele proportions in each subpopulation (p 1 , p 2 ) were randomized for each subpopulation in each replicate of each scenario (0 ≤ p ≤ 1). Each generation included stochastic binomial sampling of the available parent alleles to establish the allele proportions for the offspring, followed by deterministic symmetrical dispersal to create the parent populations for the next generation. For each scenario (combination of m, N), there were 1000 independent replicate iterations, which could be regarded either as 1000 different pairs of populations, or one pair of populations with 1000 SNP loci showing independent segregation (i.e., in linkage equilibrium); such data are now commonplace. Note that researchers using SNP data typically search for linkage disequilibrium between pairs of loci, and remove one of each pair. If researchers obtain less than 1000 SNP loci that are unlinked (Waples et al., 2022), this will reduce the precision of any differentiation measure, including G ST or A A. Note that if there is moderate to high recombination, the number of chromosomes does not directly limit the number of statistically unlinked loci. Each iteration was run for 200 generations. Because the calculations in Appendix 1 are for equilibrium given values of m, μ, and N, it was important to ensure that this number of generations was long enough to allow a close approach to equilibrium. This was ensured in three ways. First, results were inspected to ensure that each scenario had asymptoted to a stable value for BCAFD, well before the final generation. Second, iterations were each also inspected to ensure that the variance of BCAFD between-generations was much lower than variance between replicate iterations (typically one tenth or less).
Finally, 200 generations was much greater than the expected time for F ST to reach half drift-dispersal equilibrium (t 1∕2 eq generations), which for diploids is (Crow & Aoki, 1984;Whitlock, 1992): and for haploids is where symbols are as in Appendix 1. Maximum time to halfequilibrium is 34 generations for the scenarios trialed in the main paper (Table A2.1). Given that BCAFD is a function of F ST or G ST , is seems reasonable to assume that this will also approximate the time to half-equilibrium for BCAFD and A A. The simulations should be run for several times this t 1∕2 eq . For all simulated scenarios, a time of 200 generations was chosen, which is well in excess of the expected times to half-equilibrium in Table A2.1. There is a second, opposing, constraint on the number of generations. As well as the need to ensure close approach to equilibrium, the calculations in Appendix 1 assume no fixation (i.e., loss of all alleles except one), so that it was important to run the simulations for (A2.1b) t 1∕2 eq = ln0.5 ln (1 − m) 2 1 − 1 N times that are short enough to avoid fixation. This is also important because most researchers, or the companies that do their genotyping, will filter out invariant (fixed) SNPs from the data. Table A2.2 shows that it is possible to choose simulation generation numbers that are short enough to give minimum fixation, but sufficiently large to give approximate equilibrium (Table A2.1). In the case of two equal-sized subpopulations making up a metapopulation with dispersal, N for metapopulation ≈2 × N-subpopulation; for haploid we use 4 N(metapop) instead 8 N in the equation of expected time to fixation, and then we find that t fix generations is given by ( [Kimura & Ohta, 1969, Maruyama, 1970, Crow & Kimura, 1970  where symbols are as in Appendix 1. Times to fixation in generations, for the scenarios trialed in main paper, are shown in Table A2.2. This table shows that even for the smallest effective population size of 1000, fixation in about 200 generations would require a mean initial p of 0.01 or less, for example, initial p 1 = p 2 = 0.01. With random assignment of initial p values in each of the two locations, such a situation would arise in only 0.01 × 0.01 × 100% = 0.01% of replicates. In an extreme case where N for the metapopulation was equal to the N for either subpopulation, the fixation times would be halved, and yet most of the fixation times would still be orders of magnitude larger than the 200 generations simulated. In case any fixation did occur, the program included a trap for fixation, and it was designed so that if fixation occurred in any iteration, then the iteration would be replaced by restarting from generation zero, in line with the filtering normally applied to such data. Because of the relatively short number of generations (200), there were virtually no restarts for fixation.
The simulations used a binomial mechanism for transmission of alleles between generations, because of the initial focus on 2-allele SNPs. However other mechanisms such as Poisson or negative binomial might give different results (Warton & Hui, 2017), and this might be appropriate in other cases outside the scope of this paper, including where the underlying biological process for transmitting variants is different, or not adequately understood at present.
The accuracy of the simulations was tested in two ways: • By limited comparison with results of EASYPOP (Balloux, 2001), in which the starting proportions could be set to be approximately p 1 = p 2 = 0.5. These showed almost identical results for G ST when the MATLAB simulation was run with initial p 1 = p 2 = 0.5 exactly.
• By checking G ST results from the simulation against the well-   (Takahata, 1983), using 4 N in place of 8 N for the haploid loci simulated. The intercept is very close to zero as expected. The slope coefficient (0.97) has 95% confidence limits 0.95 to 1.00. The red line is the expected 1:1 relationship, indistinguishable from the actual regression line. Note that for each scenario, the coefficient of variation for simulated G ST was always larger than the coefficient of variation for A A (sometimes double).